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ABSTRACT 

Observations indicate that intermediate mass stars, binary stars, and stehar rem- 
nants often host planets; a complete explanation of these systems requires an under- 
standing of how planetary orbits evolve as their central stars lose mass. Motivated by 
these dynamical systems, this paper generalizes in two directions previous studies of 
orbital evolution in planetary systems with stellar mass loss: [1] Many previous treat- 
ments focus on constant mass loss rates and much of this work is carried out numerically. 
Here we study a class of single planet systems where the stellar mass loss rate is time 
dependent. The mass loss rate can be increasing or decreasing, but the stellar mass 
always decreases monotonically. For this class of models, we develop analytic approx- 
imations to specify the final orbital elements for planets that remain bound after the 
epoch of mass loss, and find the conditions required for the planets to become unbound. 
We also show that for some mass loss functions, planets become unbound only in the 
asymptotic limit where the stellar mass vanishes. [2] We consider the chaotic evolution 
for two planet systems with stellar mass loss. Here we focus on a model consisting 
of analogs of Jupiter, Saturn, and the Sun. By monitoring the divergence of initially 
similar trajectories through time, we calculate the Lyapunov exponents of the system. 
This analog solar system is chaotic in the absence of mass loss with Lyapunov time 
Tiy ~ 5 — 10 Myr; we find that the Lyapunov time decreases with increasing stellar mass 
loss rate, with a nearly linear relationship between the two time scales. Taken together, 
the results of this paper help provide an explanation for a wide range of dynamical 
evolution that occurs in solar systems with stellar mass loss. 

Subject headings: planets and satellites: dynamical evolution and stability — planet- 
star interactions — stars: evolution — stars: mass loss — white dwarfs 



1. Introduction 



Solar systems orbiting other stars display a diverse set of architectures and motivate further 
studies concerning the dynamics of planetary systems. Part of the richness of this dynamical 
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problem arises from the intrinsic complexity of N-body systems, even in the absence of additional 
forces (Murray & Dermott 1999). The ledger of physical behavior experienced by such systems is 
enormous, and includes mean motion resonances, secular interactions, and sensitive dependence on 
the initial conditions (chaos). Additional complications arise from additional forces that are often 
present: During early stages of evolution, circumstellar disks provide torques that influence orbital 
elements, and turbulent fluctuations act on young planets. Over longer time scales, solar systems 
are affected by tidal forces from both stars and planets, and by general relativistic corrections that 
lead to orbital precession. Another classic problem in solar system dynamics concerns planetary 
orbits around central stars that are losing mass (Gylden 1884, Jeans 1924; see also Hadjidemetriou 
1963, 1966). Although this issue has received some recent attention (see below), this paper expands 
upon existing work in two main directions. Recent work often focuses on the particular case of 
constant mass loss rates, although stellar mass loss rates typically vary with time; in addition, 
this recent work is primarily carried out numerically (note that Veras et al. 2013 use numerical 
simulations to consider more realistic, time-dependent stellar mass loss). In this paper, for single 
planet systems, we extend existing calculations to account for time dependence of the mass loss 
rates and obtain a number of analytic results. For systems with two or more planets, we also show 
that the Lyapunov exponents, which determine the time scales for chaotic behavior, depend on the 
time scales for mass loss. As outlined below, these two results can account for a great deal of the 
possible behavior in solar systems where the central star loses mass. 

A number of previous studies have considered planetary dynamics for host stars that are losing 
mass. For our own Solar System, long term integrations have been carried out to study the fate 
of the planets in light of mass loss from the dying Sun (Duncan & Lissauer 1998). Recent related 
work estimates an effective outer boundary to the Solar System (due to stellar mass loss) in 
the range = 10^ — 10'^ AU, where orbiting bodies inside this scale remain safely bound (Veras 
k, Wyatt 2012). Planets orbiting more massive stars, which lose a larger percentage of their mass, 
have their survival threatened by possible engulfment during the planetary nebula phase (Villaver 
&; Livio 2007, Mustill & Villaver 2012), and are more likely to become unbound due to stellar mass 
loss alone (Villaver & Livio 2009). In the future of our own system, Earth is likely to be engulfed 
by the Sun (Schroder &; Connon Smith 2008), but planets in wider orbits are expected to survive. 
However, gaseous planets that escape engulfment are still subject to evaporation and can experience 
significant mass loss (Bear &; Soker 2011, Spiegel &; Madhusudhan 2012). For planets orbiting stars 
that are losing mass, a more general treatment of the dynamics has been carried out for both single 
planet systems (Veras et al. 2011) and multiple planet systems (Veras & Tout 2012; Voyatzis et 
al. 2013); these studies provide a comprehensive analysis for the particular case of constant mass 
loss rates. In addition to causing planets to become unbound, stellar mass loss can drive orbital 
evolution that leads to unstable planetary systems surrounding the remnant white dwarfs remaining 
at the end of stellar evolution (Debes k, Sigurdsson 2002). Indeed, observations indicate that white 
dwarfs can anchor both circumstellar disks (Melis et al. 2010) and planetary systems (Zuckerman 
et al. 2010); many white dwarf atmospheres contain an excess of heavy elements (Melis et al. 
2010; Jura 2003), which is assumed to be a signature of accretion of a secondary body (a planet or 
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asteroid). Finally, mass loss in binary star systems can lead to orbital instability, allowing planets 
to change their host star (Krattcr & Perets 2012; for additional related work, see also Perets & 
Kratter 2012, Moeckel k Veras 2012). 

This paper builds upon the results outlined above. Most previous studies have focused on stellar 
mass loss rates that are constant in time, and most recent work has been carried out numerically. 
However, stars generally have multiple epochs of mass loss, the corresponding rates are not constant, 
and these solar systems span an enormous range of parameter space. It is thus useful to obtain 
general analytic results that apply to a wide class of mass loss functions. The first goal of this 
paper is to study single planet systems where the stellar mass loss rate varies with time. As the 
system loses mass, the semimajor axis of the orbit grows, and the planet becomes unbound for 
critical values of the stellar mass fraction rrif = Mf/Mo^,. For systems that become unbound, we 
find the critical mass fraction m/ as a function of the mass loss rate and the form of the mass loss 
function. In other systems with mass loss, the orbit grows but does not become unbound. In these 
cases, we find the orbital elements at the end of the mass loss epoch, again as a function of the 
mass loss rate and the form of the mass loss function. For initially circular orbits and slow mass 
loss (time scales much longer than the initial orbital period), the critical mass fraction and/or the 
final orbital elements are simple functions of parameters that describe the mass loss rate. For initial 
orbits with nonzero eccentricity, however, the outcomes depend on the initial orbital phase. In this 
latter case, the allowed values of the critical mass fraction ruf (or the final orbital elements) take 
a range of values, which we estimate herein. 

Next we consider the effects of additional planets on the results described above. If the planets 
are widely spaced, they evolve much like individual single planet systems. However, if the planets 
are sufficiently close together so that planet-planet interactions are important, the systems are 
generally chaotic. As a secondary goal, this paper estimates the Lyapunov exponents for two- 
planet systems with stellar mass loss. For the sake of dcfiniteness, we focus on planetary systems 
containing analogs of the Sun, Jupiter, and Saturn, i.e., bodies with the same masses and (usually) 
the same orbital elements. We find that the time scale for chaos (the inverse of the Lyapunov 
exponent) is proportional to the mass loss time scale. As a result, by the time the star has 
lost enough mass for the planets to become unbound, the planets have begun to erase their initial 
conditions through chaos. For systems that evolve far enough in time, one can use the semi-analytic 
results derived for single planet systems with initially circular orbits (see above) as a rough estimate 
of the conditions (e.g., the final value ^/ of the radius) required for a planet to become unbound. 
Multiple planets and nonzero initial eccentricity act to create a distribution of possible values (e.g., 
for ^f) centered on these results. Since the two-planet systems are chaotic, and display sensitive 
dependence on initial conditions, one cannot unambiguously predict the value of required for an 
planet to become unbound. 

For completeness, we note that the problem of planetary orbits with stellar mass loss is anal- 
ogous to the problem of planetary orbits with time variations in the gravitational constant G. For 
single planet systems (the pure two-body problem), the gravitational force depends only on the 
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(single) product GM^, so that the two problems are equivalent (e.g., Vinti 1974). However, for the 
case with time varying gravitational constant, the product GM^, could increase with time. Current 
experimental limits show that possible variations occur on time scales much longer than the cur- 
rent age of the universe (see the review of Uzan 2003), so that these effects only (possibly) become 
important in the future universe (Adams &: Laughlin 1997). We also note that when considering 
time variations of the constants, one should use only dimensionless quantities, in this case the 
gravitational fine structure constant ac = Gm?p/ch (e.g., Duff et al. 2002). 

This paper is organized as follows. We first present a general formulation of the problem of 
planetary orbits with stellar mass loss in Section [2] and then specialize to a class of models where 
the mass loss has a specific form (that given by equations [13] and [H]). These models include a 
wide range of behavior for the time dependence of the mass loss rates, including constant mass loss 
rates, exponential mass loss, and mass loss rates that decrease quickly with time; these results are 
described in Section [3l In the following Section U we consider two planet systems and calculate 
the Lyapunov time scales for a representative sample of mass loss functions. Next we apply these 
results to representative astronomical systems in Section [5l The paper concludes, in Section [H with 
a summary of our results and a discussion of their implications. 



In this section we develop model equations for solar systems where the central star loses mass. 
We assume that mass loss takes place isotropically, so that the rotational symmetry of the system 
is preserved and the total angular momentum is a constant of motion. This constraint is explicitly 
satisfied in the analytical solutions that follow. For the numerical solutions, this property is used 
as a consistency check on the numerical scheme. On the other hand, the total energy of the system 
is not conserved because the total mass decreases with time (equivalently, the system no longer 
exhibits time reversal symmetry). 

General forms for the equations of motion with variable stellar mass are presented in many 
previous papers (from Jeans 1924 to Veras et al. 2011). Some of the subtleties of the various 
approaches are outlined in Hadjidemetriou (1963). In this section and the next we specialize to 
systems with a single planet and focus on the case where the planet mass is much smaller than the 
stellar mass, Mp <C M^,. The specific angular momentum J can be written in the form 



where ao is the starting semimajor axis and Mq* is the initial mass of the star. Equation ([T|) can 
be taken as the definition of the angular momentum parameter ry. For a starting circular orbit r/ 
= 1, whereas eccentric orbits have r/ = 1 — < 1, where e is the initial orbital eccentricity. The 
radial equation of motion can be written in the dimensionless form 



2. Model Equations for Orbits with Stellar Mass Loss 



j2 = GMo^aov 



(1) 



T] m{t) 



(2) 
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where r] is constant and where we have defined a dimensionless (total) mass 

The dimensionless radial coordinate ^ = r/ao and the dimensionless time variable is given in units 
of 17-1 ^ (a3/GMo*)i/2. 

The goal of this work is to find general solutions to the problem where the dimensionless mass 
m{t) monotonically decreases with time. In the reduction of the standard two-body problem to an 
analogous one-body problem, the equation of motion describes the orbit of the reduced mass. In 
the version of the problem with mass loss represented by equation ([2]), the motion is also that of a 
reduced mass (e.g., Jeans 1924, MacMillan 1925, Hadjidemetriou 1963). In this treatment, the mass 
loss functions (defined below) refer to the dimensionless scaled mass m{t). The model equation ([2]) 
is exact in the limit where the planet mass is small compared to the stellar mass, i.e., Mp/M^ — t- 0. 
For finite planetary masses, the scaling between the two-body problem and the equivalent single 
body problem changes quantities by 0{Mp/M^). In applications of interest, however, our choice of 
mass loss functions (and their uncertainties) provides the greatest degree of approximation - much 
larger than than the 0{Mp /M^) difference between the reduced problem and the full problem. 

This paper thus focuses on the dimensionless problem of equation ([2]). The starting semimajor 
axis is unity (by definition) and the initial conditions require that the starting radial coordinate 
lies in the range 1 — e<.^o ^ 1 + e, where the starting eccentricity e is given hy rj = 1 — e^. Note 
that choosing the value of is equivalent to specifying the starting phase of the orbit (up to a 
sign). The initial energy £q = —1/2 by definition and the initial (radial) velocity .^o is given by 

/2 _ 2Co-r?-Cg _ [(l + e)-go][Co-(l-e)] 

^0 - 72 - 72 • 

so so 

The initial velocity can be positive or negative, where the choice of sign completes the specification 
of the starting phase of the orbit. 



2.1. Change of Variables 

The equation of motion ([2]) is complicated because it contains an arbitrary function, namely 
m{t), that describes the mass loss history. On the other hand, the independent variable (time) 
does not appear explicitly. As a result, we may define a new effective "time" variable u through 
the expression 

u=-, (5) 
m 

where m = m{t). The generalized time variable u starts at n = 1 and is monotonically increasing. 
In terms of the variable u, the basic equation of motion ^ takes the equivalent form 

od^C dS, r] 1 



- 6 - 



Next we note that both standard lore and numerical solutions (beginning with Jeans 1924) 
show that, in physical units, the product aM « constant. In terms of the current dimensionless 
variables, this finding implies that the function 

/ = - = em (7) 

u 

should vary over a limited range. We thus change the dependent variable from ^ to / and write 
the equation of motion in the form 



lnY' + 2n/') + ^(n/' + /) 



r]-f, (8) 



where primes denote derivatives with respect to the variable u. Keep in mind that equation ([8]) 
is equivalent to the original equation of motion with a change in both the independent and 
dependent variables. 

The leading coefficient in equation ([8]) represents an important quantity in the problem: Note 
that the time scale for mass loss is given hy u/ii and the orbital time scale is given by v? f^^"^ (this 
latter time scale is the inverse of the orbital frequency, and is shorter than the orbital period by a 
factor of 2-k). The ratio A of these two fundamental time scales is given by 

A2 ^ uWf . (9) 

The leading coefficient in equation ([8]) is thus A^, the square of the ratio of the orbital time scale 
to the time scale for mass loss. For small values of A^, the mass loss time is long compared to 
the orbit time, and the orbits are expected to be nearly Keplerian; for larger A'^, the star loses a 
significant amount of mass per orbit and a Keplerian description is no longer valid. For the former 
case, where mass loss is slow compared to the orbit time, we can use the parameter A to order the 
terms in our analytic estimates. 

In addition to the coefficient A^, given by the ratio of time scales, another important feature 
of equation ([8]) is the index /3 appearing within the square brackets, where 

(10) 

The index (5 encapsulates the time dependence of the mass loss. This paper will focus on model 
equations with constant (3 (such models have a long history, from Jeans 1924 to Section [2.2p . 

The Orbital Energy: For the chosen set of dimensionless variables, the energy £ of the system 
takes the form 

The energy has a starting value £ = —1/2, by definition, and increases as mass loss proceeds. If 
and when the energy becomes positive, the planet is unbound. Although the energy expression 
(fTTI) appears somewhat complicated, the time dependence of the energy reduces to the simple form 

d£ _ 1 
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Note that the derivative of the energy is positive definite, so that the energy always increases. 
Since the energy is negative and strictly increasing, the semimajor axis of the orbit, when defined 
according to a oc |<?|~^, is also monotonically increasing. 



2.2. Mass Loss Functions 

Next we want to specialize to the class of mass loss functions where /3 = constant. The defining 
equation (|lUp for the mass loss index can be integrated to obtain the form 

u = juf^, (13) 

where 7 is a constant that defines the mass loss rate at the beginning of the epoch (when t = 0, 
m = 1, and u = 1). For a given (constant) value of the index /3, the dimensionless mass loss rate 
has the form 

m = -^rn^^-P) . (14) 

In addition to simplifying the equation of motion, this form for the mass loss function is motivated by 
stellar behavior, as discussed below. The dimensionless parameter 7 is defined to be the ratio of the 
initial orbital time scale to the initial mass loss time scale. Specifically, if we define r = (M^,/M*)o, 
then 7 is given by 

1/2 

(15) 

1.6x10-/ - W "0 
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where Mq* is the stellar mass and oq is the semimajor axis at t = 0. For typical orbits, ao = 1 - 
100 AU, so that we expect the parameter 7 to be small, often in the range 10^ ^ 7 ^ 10~^. 

The mass loss rate of stars is often characterized by the physically motivated form 

where Mq is constant and depends on the phase of stellar evolution under consideration (Kudritzki 
&: Reimers 1978, Hurley et al. 2000). Since the radius and luminosity depend on stellar mass 
(for a given metallicity), the physically motivated expression of equation (jl6p can take the same 
power-law form as equation (fHl) . which corresponds to a constant mass loss index (see equations 
[ID] and [S]). 

Using the scaling law ()16p . the power-law index appearing in equation ()14p can be positive or 
negative, depending on how the stellar luminosity and radius vary with mass during the different 
phases of mass loss (see Hurley et al. 2000 for a detailed discussion). For example, if we consider 
main-sequence stars, the stellar cores adjust quickly enough that the luminsoity obeys the standard 
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mass-luminosity relationship L^, ~ M* (where the index p ~ 3) and the mass-radius relationship 
i?* ~ MJ' (where the index q typically falls in the range 1/2 < q < 1). For main-sequence stars we 
thus obtain the scaling law m ~ — m"™ , where the index am = p + q — 1 is predicted to lie in the 
range 2.5 < am < 3; the corresponding mass loss index lies in the range —1 < /3 < —1/2. Next 
we consider stars on the first giant branch or the asymptotic giant branch. In this phase of stellar 
evolution, mass loss occurs from an extended stellar envelope, but the luminosity is produced deep 
within the stellar core. As the star loses mass, the core and hence the luminosity remains relatively 

— 1/3 

constant, whereas the radius scales approximately as i?* ~ ' (Hurley et al. 2000). For this 
case, one obtains the scaling law m ~ —m~^/^, with a mass loss index /3 = 10/3. In general, for 
m ~ — m"™, the mass loss index (3 = 2 — am- As these examples show, the mass loss index can 
take on a wide range of values — 1 < /? < 4. 

To fix ideas, we consider the time dependence for mass loss functions that are often used. For 
a constant mass loss rate, the most common assumption in the literature, the index (3 = 2, and the 
mass evolution function has the form 

m{t) = 1 - 7t and u{t) = ■ (1"^) 

The value (3 = 2 marks the boundary between models where the mass loss rate accelerates with 
time (/3 > 2) and those that decelerate {(3 < 2). For the case of exponential time dependence of 
the stellar mass, the index (3 = 1, and the mass loss function has the form 

m(t) = exp[— 7t] and u(t) = exp[7t] . (18) 

The value (3 = 1 marks the boundary between models where the system reaches zero stellar mass 
in a finite time {(3 > 1), and those for which the mass m — )• only in the limit u — )• cx). For the case 
with index (3 = 0, which represents an important test case, the mass evolution function becomes 

m{t) = Y^~^ ^^'^ = 1 + 7t . (19) 

For (3 = 0, analytic solutions are available (see Section [3. 2p . which inform approximate treatments 
for more general values of the index (3. Finally, the case where (3 = —1 plays a defining role (see 
Section 13. 1|) and corresponds to the forms 

m{t) = (1 + 27t)-i/2 and u{t) = (1 + 27*)^/^ . (20) 

The value (3 = —1 marks the boundary between models where the planet becomes unbound at finite 
stellar mass (/3 > — 1) and those for which the planet becomes unbound only in the limit m — t- or 
u — 7- 00 {(3 < —1). In general, for constant /3 7^ 1, the time dependence of the mass takes the form 

m{t) = -^^=[l-{(3-lht]'/^^~'^ . (21) 

The particular case (3 = 1 results in the decaying exponential law of equation ()18p . 
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2.3. Equation of Motion with Constant Index (3 

For constant values of the mass loss index /3, the equation of motion reduces to the form 

a2 [u^f" + (2 + p)uf' + I3f]=r]-f, (22) 
where the ratio of time scales A is given by 

A2 = ^2^2/3+2 j3 (23) 

By writing the equation of motion in the form ()22p . we immediately see several key features of the 
solutions: 

When the parameter A <C 1, the left-hand side of equation (j22p is negligible, and the equation 
of motion reduces to the approximate form / « ry = constant. This equality is only approximate, 
because the function / also executes small oscillations about its mean value as the orbit traces 
through its nearly elliptical path (see below). Nonetheless, this behavior is often seen in numerical 
studies of planetary systems with stellar mass loss (e.g., see Debes & Sigurdsson 2002). Orbital 
evolution with A <C 1 is often called the "adiabatic regime". We note that this terminology is 
misleading, however, because "adiabatic" refers to evolution of a thermodynamic system at constant 
energy (heat), whereas the systems in question steadily gain energy through stellar mass loss (the 
gravitational potential becomes less negative). 

When the parameter A ^ 1, the left-hand side of equation ()22p dominates, and the solutions 
for f{u) take the form of power-laws with negative indices. In this regime, the equation of motion 
approaches the form 

nY' + (2 + /3)^z/' + /3/ = 0, (24) 
so that the function f{u) has power-law solutions with indices p given by the quadratic equation 

{p+l){p + p) = 0. (25) 

The general form for the solution f{u) in this regime is thus 

fin) = - + 4 • (26) 

After the solutions enter this power-law regime, the energy can quickly grow and the planet can 
become unbound. To illustrate this behavior, consider the differential equation ()12p for the orbital 
energy. We first consider the regime where A <C 1 and the function / is nearly constant. For the 
benchmark case / = 1, the equation can be integrated to obtain 

As long as / ~ 1, the energy remains negative and the planet remains bound, except in the limit 
ti — 7- oo. Now let A > 1 so that the solutions enter into the power-law regime. If we let the solution 
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have the form / = Aju, for u > Uc, the differential equation for the energy can be integrated to 
obtain 

where the subscript c denotes the reference point where the solutions enters into the power-law 
regime. Since Auc ~ u'^ and \£c\ ~ '"c/^j energy quickly becomes positive once the power-law 
regime is reached. This argument indicates that the planet becomes unbound when the time scale 
ratio A is of order unity (as seen in previous studies) . The subsequent subsections provide further 
verification of this finding. 

In order for the solution to make the transition from / ~ constant to the power-law solutions 
that cause the orbits to become unbound, the ratio of time scales A must grow with time. However, 
growth requires that /? > — 1 (see equation j23j). We can understand this requirement as follows: 
The orbital time scale P varies with u (and hence time) according to P ~ u^f^^'^. Since / is nearly 
constant, this relation simplifies to the form P ^ . The time scale for mass loss r is given by 
r = u/ii, which has the form r ~ v}~^ from equation (jl3|) . As a result, when /3 = — 1, the orbit 
time has the same dependence on stellar mass as the mass loss time scale, so that the ratio A is 
nearly constant as the star loses mass. For (5 < —1, the ratio A of time scales decreases with time, 
and the system grows "more stable". 



3. Results for Single Planet Systems with Stellar Mass Loss 

This section presents the main results of this paper for single planet systems with a central 
star that loses mass. First, we consider mass loss index /3 = — 1, which marks the critical value such 
that systems with /3 > — 1 become unbound at finite values of the stellar mass, whereas systems 
with (3 < —1 only become unbound in the limit m — )• 0. Next we consider mass loss index /3 = 0; in 
this case, the solutions can be found analytically, and these results guide an approximate analytic 
treatment of the general case, which is addressed next. We also consider the limiting case where 
stellar mass loss takes place rapidly. 



3.1. The Transition Case 

Here we consider systems where the mass loss index /3 = — 1, which corresponds to the tran- 
sition value between cases where the ratio A of time scales grows with time (/3 > —1) and those 
where the ratio decreases with time (/3 < — 1). In this regime, the equation of motion reduces to 
the form 

7^[.Y' + n/'-/]=-^--i. (29) 
The equation of motion can be simplified further by making the change of (independent) variable 



w = log u , 



(30) 
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so that the equation of motion becomes 



7 



f 



P 



1 

P 



(31) 



This version of the equation of motion (first considered by Jeans 1924) contains no exphcit de- 
pendence on the independent variable w, so that the equation can be integrated to obtain the 
expression 

where is a constant that plays the role of energy. Note that we have chosen the sign such that 
E > {) and that E = 1 for initially circular orbits. In order for the function f{w) to have oscillatory 
solutions, the fourth order polynomial 



Pif) 



Ep + 2f-rj 



(33) 



must be positive between two positive values of /. In order for this requirement to be met, the 
parameters must satisfy the inequalities 

9 



rjE < 



and 



^<{E/3f/'. 



(34) 



The first inequality is always satisfied for the cases of interest. The second inequality in equation 
P4p determines the maximum value of the mass loss parameter 7 for which oscillatory solutions 
occur. 



3.2. Systems with Vanishing Mass Loss Index 

In the particular case where /? = 0, the equation of motion can be simplified. In particular, 
the first integral can be taken analytically to obtain the form 

7^(-V)' = -;^ + |-i?, (35) 

where E is a constant of integration. The parameter E plays the role of energy for the orbit problem 
where the function f{u) plays the role of the radial coordinate. Although E is constant, the energy 
S of the physical orbit (where is the radial coordinate) increases with time. Notice also that we 
have adopted a sign convention so that E > 0. The value of E depends on the initial configuration. 
For the particular case where the orbit starts at periastron, for example, the initial speed ^ = 
and the energy constant has the value 

E=l-P{l-ef, (36) 

where the eccentricity e = (1 — r])^^"^. In general, the initial value /o = ^^o can lie anywhere in the 
range 1 — e</o<l + e, and the energy constant has the general form 

E = l- Pp, ± 27 [2/0 - 7? - /2] , (37) 
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where the choice of sign is determined by whether the planet is initially moving outward (+) or 
inward (— ). With the energy constant E specified, the turning points for the function / are found 
to be 

If we consider the function / to play the role of the radial coordinate, then equation (j38p defines 
analogs of the semimajor axis a^, and eccentricity e^,, which are given by 

a* = — and = [1 — r]E] ' . (39) 

E 

For a given starting value /o, the effective eccentricity is given by the expression 

el = e'+ 7^(1 - e')f, T 2j{l - e') (2/o - v " fo) ■ (40) 

Note that the effective eccentricity of the function / is larger than the initial eccentricity e of 
the original orbit (before the epoch of mass loss). In particular, for a starting circular orbit e = 0, 
the effective eccentricity e* = 7 7^ 0. The integrated equation of motion ([55]) can be separated and 
written in the form 

fdf ^1/2 du 



[(/-/i)(/2 -/)]'/' 7 



u 



2 ■ 



(41) 



If we integrate this equation from one turning point to the other, the change in mass of the system 
the same for every cycle, i.e.. 

Am = , (42) 

where E is given by equation ()37p . Following standard procedures, we can find the solution for the 
orbit shape, which can be written in the form 

1 = !l^ = l + e^cos9. (43) 
J ? 

The orbit equation thus takes the usual form, except that the original eccentricity e is replaced 
with the effective eccentricity e* and the effective "radius" variable (/) scales with the mass/time 
variable u = 1/m. 

For this mass loss function (with /? = 0), we can find a simple relationship between the value of 
the time scale ratio A and the value of / when the planet becomes unbound. To obtain this result, 
we insert the first integral from equation (I35p into the general expression (llip for the energy £ of 
the orbit and set £ = 0. After eliminating the derivative /', we can solve for the time scale ratio 
Aj as a function of the final value of /. When the planet becomes unbound, the time scale ratio is 
thus given by 

,^^ i2f-.r-M^f-.-r-f\ ,,,, 

Here, / is evaluated when the planet becomes unbound. In this case, however, the value of / is 
constrained to lie in the range /i < / < /2, where the turning points are given by equation ()38p . 
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For small 7, the orbit oscillates back and forth between the turning points many times before the 
planet becomes unbound. The final value of / is thus an extremely sensitive function of the starting 
orbital phase. This extreme sensitivity is not due to chaos, and can be calculated if one knows the 
exact orbital phase at the start. In practice, however, the final value of / can be anywhere in the 
range /i < / < /2. 

Figure 13.21 shows the final values of the time scale ratio A as a function of the final value of 
f = i/u = ^m. Curves are shown for nine values of the starting eccentricity, where e = 0.1, 0.2, 
. . . 0.9. The innermost (outermost) closed curve in the figure corresponds to the smallest (largest) 
eccentricity. Each value of / corresponds to two possible values of the time scale ratio A, one for 
orbits that are increasing in / and one for orbits that are decreasing in / at the time when the 
planet becomes unbound. These two values of A (for a given /) correspond to the two branches of 
the solution given by equation (l44l) . 

For comparison. Figure 13.21 shows the same plane of parameters for the final values of the 
time scale ratio \f and ff for planetary systems with constant mass loss rate (where /3 = 2). 
These results were obtained through numerical integration of equation (j22p . Note that the range 
of allowed values for the time scale ratio Aj is much larger than for the case with /? = 0, whereas 
the range of final values // is somewhat smaller. 

One can also show that for circular orbits [r] = \ and e = 0) , the value of the time scale ratio 
A = 1 when the planet becomes unbound. For circular orbits in the limit 7—7-0, the turning points 
of the orbit appraoch /1.2 = 1. Using / = 1 and r/ = 1 in equation (j44p . we find A = 1. 



3.3. Limit of Rapid Mass Loss 

If we now consider the case where the mass loss is rapid, so that the equation of motion has 
solution (f26ll throughout the evolution, we can fix the constants A and B by applying the initial 
conditions. Since / = S^/u and « = 1 at the start of the epoch, /(I) = ^1, where ^1 is the starting 
value of the orbital radius. By definition, the semimajor axis is unity, and the starting orbital 
eccentricity is given by = 1 — ry. The starting radius thus lies in the range 

i-yr^<6<i + \/r^, (45) 

which is equivalent to 1 — e < ^1 < 1 + e. The derivative /' = df /du is given by 

/' = -^ + ^ = -^ + ^^. (46) 
V? u v? 7ii^+i dt 

In the regime of interest where 7^1, the second term is small compared to the first. In this limit, 
/'(I) = —^1, where .^1 lies in the range indicated by equation ([l5|). The constants A and B are 
thus determined to be B = and ^ = ^1, so that the solution has the simple form 



f{u) = ^ . (47) 
u 
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Fig. 1. — Time scale ratio A/ as a function of the final value of / = // for planetary systems with fi 
= 0. For a given angular momentum ry, specified by the starting eccentricity, the allowed values of 
\f form closed curves in the plane as shown. Curves are shown for a range of starting eccentricity, 
from e = 0.1 (inner curve) to e = 0.9 (outer curve). 
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Fig. 2. — Time scale ratio A/ as a function of the final value of / = // for planetary systems 
with (3 = 2 (constant mass loss rates). For a given angular momentum t/, specified by the starting 
eccentricity, the allowed values of Aj form closed curves in the plane as shown. Curves are shown 
for a range of starting eccentricity, from e = 0.1 (inner curve) to e = 0.9 (outer curve). Compare 
with Figure 13.21 (and note the change of scale) . 
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The energy of the orbit, by definition, starts at £{u = 1) = £i = —1/2, and the energy obeys the 
differential equation (jl2p . Combining the solution of equation ()47p with the differential equation 
()12p for energy, we can integrate to find the energy as a function of u (equivalent to time or mass) , 

We can then read off the value of Uf, and hence the mass mj, where the energy becomes positive 
and the planet becomes unbound, i.e., 

», = i- = l-|. (49) 

Note that this critical value of the mass depends on the orbital phase of the planet within its orbit, 
i.e., the result depends on rather than the starting semimajor axis, which is unity (notice also 
that this condition is equivalent to that given by equations [46-48] in Veras et al. 2011). For circular 
orbits, we must have = 1, so that planets become unbound when the stellar mass decreases by 
one half (as expected). 

For cases where the mass loss is rapid, but the planet remains bound, we can find the orbital 
properties for the post-mass-loss system. Consider the limiting case where the star has initial mass 
m = 1, and loses a fraction of its mass instantly so that it has a final mass moo, i-e., 

m{t) = moo + (1 - moo)H{-t) , (50) 

where H is the Heaviside step function. The mass loss thus occurs instantaneously at t = 0. For 
t < 0, the solutions to the orbit equation ([2|) have the usual form. 



,51) 



where .^2,1 = 1 ± e and rj = 1 — e^. After mass loss has taken place, the new (dimensionless) stellar 
mass is moo, and the first integral of the equation of motion can be written in the form 



■2 _ 2moo rj ^ 2(1 - mpp; 

, (52) 



where the final constant term takes into account the change in (dimensionless) energy at the moment 
of mass loss. The radial position at t = is ^oi since the planet is initially in a bound elliptical 
orbit, the radial coordinate must lie in the range 1 — e<^o^l + e. The energy £f of the new 
orbit is thus given by 

2(1 — m„o) / X 

2£f = -1 + ^ ^ ' . (53) 

The energy Sj is negative, and the orbit is bound, provided that the remaining mass moo > 1 — Co/2. 
This condition is thus consistent with equation (jM]), which defines the mass scale at which orbits 
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become unbound in the limit of rapid mass loss. In terms of the energy iS'j, the turning points of 
the new orbit take the form 



i± = ^ ' ^j^^i ' "'^ . (54) 
We can then read off the orbital elements for the new (post-mass-loss) orbit, i.e., 



«/ = and e/ = ^Jl - 2\£f\rj/ml^ , (55) 



where the new orbital energy £f is given by equation 



3.4. Systems with General Mass Loss Indices 

In order to address the general case, we first change variables according to the ansatz 

X = u" where a = (3 + 1 . (56) 
After substitution, the equation of motion becomes 

j^a^x^ [xVxx + 2xf^] + 72/3x2/ = ^-L, (57) 

where the subscripts denote derivatives with respect to the new variable x. The ratio of time scales 
is now given by 

A2 = -f^x^f . (58) 
We can integrate the differential equation (j57p to obtain the implicit form 



T^a^ [x^f^] 2 + 27^/3 x'^ff^dx = -^ + l-E, (59) 
where E > has the same meaning as before. To move forward, we define the integral quantity 

J = 2^^I3 r x^fUdx, (60) 



Ji 

so that 

^W[x'U]' + J=^-l^^L^. (61) 

Note that J = 0{X'^), which means that J will be negligible for most of the mass loss epoch (see 
Appendix [A|) . The energy of the system can be written in the form 



u'£ = \i^x\axh + ff + Jp--f (62) 



At the start of the evolution (where u = 1 and x = 1), the energy £ = —1/2 by definition. Using 
this specification, we can find the value of the integration constant which takes the form 

= 1 - 7^/2 ± 27 (2/o - r? - /o') ^'^ , (63) 
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where /o(= Co) is the starting value of the function (radial variable). 

At an arbitrary time during the epoch of mass loss, we can write the derivative fx = df /dx in 
the form 

7a [x'U] =±j{2f-^- Ef - Jff' . (64) 

Since J = 0{X'^), the condition \J\ <^ E holds for most times. As a result, working to leading 
order, we can set J = in equation (|64p and recover analogs to the orbital solutions found earlier 
in Section [312] (for a related result, see Radzievskii &: Gel'Fgat 1957; for an alternate approach, see 
Rahoma et al. 2009). The only difference is that the dependent (time-like) variable u is replaced 
with X = u^~^^. As a result, the turning points for the function f(x) will be given by equation (j38p 
and the orbital elements for /(x) are given by equation (p9|) . 

The basic behavior of the orbit is illustrated by Figure [331 The function /, plotted here versus 
u as the solid black curve, oscillates between the turning points (marked by the red horizontal 
lines) given by equation p8|) . The radial coordinate (here log ^ is plotted as the dotted blue curve) 
oscillates also, but grows steadily. The eccentricity of the orbit (green dashed curve) also oscillates, 
but grows with time. Finally, the time scale ratio A (magenta dot-dashed curve) also oscillates and 
grows with time. The simple oscillatory behavior for f{u) ceases near the point where the time 
scale ratio A becomes of order unity. Note that the function / falls outside the boundary marked 
by the turning points near u = 12 in the Figure. 

Next we consider the time evolution of the energy of the orbit. After some rearrangement, the 
energy from equation (j62p can be rewritten in the form 

2u'^£ = -E- J (65) 

±27x (2/ - ry - Ef - J ff'^ + ^i^x^f . 
Since J = 0{X'^), it is often convenient to work in the limit J — )• where the energy becomes 

2u^£ = -E± 2jx (2/ - ry - Ef) + j'^x'^f + O (J) . (66) 

This form for the energy shows why the product am of the effective semimajor axis and the mass 
is slowly varying: In dimensionless units, the energy £ = —m/2a so that 

(am)"^ = -2£u^ = E + O {\) . (67) 

The product am is thus nearly constant as long as the time scale ratio A is small, and the departure 
is of order A. When A is small, the orbit cycles through many turning points before the mass 
changes substantially, so that the average of the above equation becomes 

{{amY^) = {-2Eu^) = E + (A^) , (68) 

so that the average of am is constant to second order. 
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Number of Cycles: If we ignore J for now and integrate equation (|59p over one cycle, we obtain 



7a f fdf f'^ dx 1 1 



E^/^Jl (/-/l)l/2(/2-/)l/2 J, X, X2- 

The integral on the LHS gives us ir/E. If we integrate over cycles we obtain the expression 

-faNn = E^/^ [I - m%] , (70) 

where we assume that m = 1 at the start. The total number of possible cycles occurs when — t- 0, 
so that 

Nt = ^= f.^.. . (71) 
Ti^a 7r7(p + 1) 

Since £' ~ 1 and 7r(/3 + 1) ~ 10, whereas 7 <C 1, we expect the number of cycles Nt ~ 1/(107) to 
often be large. 

The Last Cycle: The above analysis (if we continue to work in the regime J ^ 1) suggests 
that the last cycle occurs when the right hand side of equation ()69p is no longer large enough to 
balance the left hand side, which is the same for each cycle. As a result, a minimum mass must be 
left in the star in order for the orbit to complete a cycle (in the function /). This condition can be 
written in the form 

which holds for a = 1 + /3 > 0. At the point when the mass falls below this threshold, the time 
scale ratio is given by A = {E f^/"^ / {tto). Since Ef ~ 1 and vra ~ 3 — 10, the system crosses the 
threshold so that / cannot complete a cycle just before the time scale ratio A reaches unity. 

Final States: If we set the energy equal to zero and replace the variable x with the time scale 
ratio Xf (evaluated at the moment that planet becomes unbound), we find the condition 

(2/ - r/)^/^ = Xff'/' ± (2/ - - Ef' - Jf^" , (73) 
which can then be written in the form 

, {2f-vy/'±{2f-r^-Ef-JfY' 

A/ = jYJ-2 • (74) 

This expression is thus a generalization of that obtained for the special case with /3 = (see equation 
|44j and Figure 13. 2p . The differences are that we have included the extra term J and that the result 
is written in terms of the variable x instead of u. 

The orbital eccentricity, calculated the usual way, oscillates with time with an increasing 
amplitude of oscillation (e.g., see Figure [3^ . As shown here, however, the function / executes 
nearly Keplerian behavior, with nearly constant turning points, where this statement is exact in 
the limit J — t- 0. The oscillation of eccentricity, although technically correct, is misleading. The 
turning points of the orbit (in the original radial variable ^) are strictly increasing functions of 
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Fig. 3. — Evolution of the orbit during the epoch of stellar mass loss. In this example, the mass loss 
function has index (3 = 2, corresponding to a constant mass loss rate. The other parameters are 7 = 
10~^, e = 0.3, and /o = 1 (initially going inward). The black curve shows the function f{u) = ^/u; 
the red horizontal lines mark the analytically determined turning points of the function. The blue 
dotted curve shows the evolution of the radial coordinate ^ (plotted here as log^oi?])- The magenta 
dot-dashed curve shows the evolution of the time scale ratio A. Finally, the green dashed curve 
shows the eccentricity e of the orbit. 
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time. The oscillation in e arises because the orbits are not ehipses, and, in part, because the period 
of the orbits in ^ are not the same as the period of the orbits in /. As a result, the osciUations 
in the calculated eccentricity do not imply that the near-elliptical shape of the orbit is varying 
between states of greater and lesser elongation. Instead, these oscillations imply that if the mass 
loss stops and the orbit once again becomes an ellipse, the value of the final eccentricity of that 
ellipse oscillates with the ending time of the mass loss epoch. 

Orbital Elements during and after the Mass Loss Epoch: Next we consider the case where the 
planet remains bound after the epoch of stellar mass loss. In this case, we want to estimate the 
orbital elements of the planet. We are mostly interested in specifying the orbital elements at the 
end of the mass loss epoch, but we can also evaluate them at any time while the star continues to 
lose mass. Suppose that the orbit passes through N turning points of the function / during the 
mass loss epoch. The orbit will then complete a partial cycle so that the final value of / lies between 
the turning points, /i < / < f2- In the ideal case, where we have complete information describing 
both the starting orbital elements and the final value of the stellar mass, and where the mass loss 
function is is described exactly by a model with constant index /?, we can calculate the final value 
ff. In practice, we will often have incomplete information: The number of cycles is generally large, 
N ^ 1 (see equation [71]), and it is unlikely that stellar mass loss can be exactly described by 
a model with constant /3 for a precise number A'' of cycles (and then stops abruptly). As result, 
we are unlikely to know where the final value of / lies between the turning points. Additional 
planets, or other perturbations, increase this uncertainty (see Section H]). In this case of incomplete 
information, we can write the mean value of the energy (averaged over the cycles) in the form 

2n^£ = -E + :^, (75) 

where we have replaced / with the value 1/E of its effective semimajor axis, and where the remaining 
term averages to zero. This estimate for the final energy has an uncertainty — a range of possible 
values — due to the lack of knowledge of where the planet lies in its orbit during the final cycle. 
This range of energy is given by the form 

AS 2jx^f-ri-Efy^ 
, 27x6* 



E3/2 _^2^2/£;3/2 ' 

where is defined by equations (|39]) and (f40l) . With the energy specified, the final value a/ of the 
semimajor axis is given by 

The expected value of the energy £f is given by equation ([75]) . but it can take on any value in the 
range defined by equation ()76p . Similarly, the final value e/ of the orbital eccentricity is given by 

e} = l + r]{2u^£f) ^ 1 - rjE + rj-f^x^/E'^ , (78) 
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where the second (approximate) equahty holds for the mean value of energy given by equation ()75p . 
Since the energy £f can have a range of values, given by equation ()76p . the orbit has a corresponding 
range of possible eccentricities. 

The expressions derived above for the final orbital elements are expected to be valid provided 
that the time scale ratio A is small compared to unity (and hence \J\ <C 1). The time evolution 
of the orbital elements is illustrated in Figure 13.41 for a representative system with mass loss index 
(3 = 2 and mass loss parameter 7 = 10~^. Numerical integration of the full equation of motion 
(solid curves) show that the semimajor axis and eccentricity both oscillate and (on average) grow 
with time. (Note that the Figure plots log[a].) The values of the elements (a, e) calculated from 
the average energy (from equation [75]) provide a good approximation to the mean evolution of 
the orbital elements (see the central dotted curves in Figure [3^ . Furthermore, using the range 
of allowed energy calculated from equation (f76]l . we can calculate upper and lower limits to the 
expected behavior of the semimajor axis and eccentricity (shown as the upper and lower dotted 
curves). Note that the solutions for a and e oscillate back and forth between these limiting curves. 
In this example, the planet becomes unbound near u = 28. Prior to that epoch, near u ~ 20, the 
limits on the energy allow for the planet to become unbound, and the upper limit for the semimajor 
axis approaches infinity. The approximation scheme thus breaks down at this point. 



3.5. Numerical Results 

The equations of motion can be numerically integrated to find the value of the radial 
coordinate when the system becomes unbound (when the energy becomes positive). For the case 
of exponential mass loss, /3 = 1, the result is shown in Figure 13.51 as a function of the mass loss 
parameter 7 (top panel). The figure shows curves for initially circular orbits (t/ = 1, smooth curve) 
and for nonzero starting eccentricity (r/ = 0.9, rapidly oscillating curve). Note that the ij = 0.9 
curve is shown only for 7 > 0.003; for smaller values of 7, the curve oscillates more quickly as a 
function of 7, and the curve would appear as a solid black band in the figure. The bottom panel 
shows the value of A, the ratio of the orbital period to the mass loss time scale, evaluated when the 
system becomes unbound. As expected, the parameter A is of order unity when the system energy 
becomes positive and the planet becomes unbound. For the case of circular orbits at the initial 
epoch (rj = 1), the value of A ~ 1.3 for small 7. For starting orbits with nonzero eccentricity, the 
value of A takes on a range of values, but remains of order unity. For the case shown [r] = 0.9, 
oscillating curve), the parameter A varies between about 1/2 and 2. 

The above trend holds over a range of values for the mass loss index /3. Figure [331 shows the 
value of the time/mass variable u = 1/m when the planet becomes unbound as a function of the 
mass loss parameter 7. Results are shown for mass loss indices in the range < /3 < 3. For all 
values of the index /3, the curves become nearly straight lines in the log- log plot for small values of 
7, which indicates nearly power-law behavior of the form ~ 7~^, where the index p = l/(/3 + 1). 
We can find a simple fitting function for the final value n/ of the inverse mass variable (for initially 




Fig. 4. — Evolution of orbital elements during the epoch of stellar mass loss. In this example, 
the mass loss function has index (3 = 2. The other parameters are 7 = 10-4, e = 0.3, and 
/o = 1 (initially going inward). The solid curves show the semimajor axis (log a; top) and orbital 
eccentricity (e; bottom), calculated from numerical integration of the equation of motion. For each 
orbital element, the three dotted curves show the average value, the upper limit, and the lower 
limit, as calculated from the analytic expressions given by equations (|75l -[78 |) . 
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Fig. 5. — Radial coordinate of the planet at the moment when the system becomes unbound, shown 
here as a function of the mass loss parameter 7 for exponential mass loss (top panel). The nearly 
monotonic curve shows the result for the case of circular starting orbits; the oscillating curve shows 
the result for angular momentum parameter r/ = 0.9 (which corresponds to starting eccentricity e 
= \/oT 0.316 . . .). Bottom panel shows the value of the parameter A when the planet becomes 
unbound, for both circular starting orbits (smooth curve) and eccentric orbits {rj = 0.9; oscillating 
curve). All orbits are started at periastron. 
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circular orbits) as a function of 7: 

, (79) 

where cq is a constant. For cq = 0.74212, the error is less than 0.4% for mass loss indices in the range 
< /9 < 4. Power-law fits resulting from equation (j79p are shown as the dashed lines in Figure [331 
The power-law fits provide a good approximation for 7 < 0.1. For larger values of 7 > 0.1, the final 
value Uf ^ 2. Equation ()79p describes the final mass values (for initially circular orbits) for realistic 
mass loss prescriptions: Most of the stellar mass is usually lost on the asymptotic giant branch 
(AGB), where the mass loss function has (3 ~ constant (see Fig. 13 of Veras et al. 2011). Note 
that mass loss on the AGB takes place through a series of pulses, but this complication primarily 
affects tides (Mustill & Villaver 2012), rather than the overall mass loss profile. 

A related result in shown in Figure 13.51 which plots the values of the ratio A of time scales, 
evaluated at the moment when the planet becomes unbound, as a function of the mass loss param- 
eter 7. In the limit of small 7, the time scale ratio A approaches a constant value (of order unity). 
The finding that A has a value of order unity (in the limit of small 7) when the planet becomes 
unbound is expected: In physical terms, this result means that the mass loss time scale has become 
shorter than the orbital period, so that the potential well provided by the star is changing fast 
enough that the orbital motion does not average it out. 

The limiting values of the time scale ratio A are shown in Figure 13.51 as a function of the mass 
loss index (3. Here the time scale ratios are evaluated at the moment when the planet becomes 
unbound. Results are shown for the limiting case of small 7 (from Figure [331 we see that the time 
scale ratio A approaches a constant value as 7 — ?• 0). All of the values are of order unity; for the 
particular case where (3 = 0, the final value of the time scale ratio A = 1. Since this function Xf{(3) 
is useful for analysis of orbits in systems losing mass, we provide a simple fit. If we choose a fitting 
function of the form 

logA^ = Cl/3 + C2/3^ (80) 

where ci and C2 are constants, we obtain a good fit with the values ci = 0.21658 and C2 = 0.04102. 
The fitting function is shown as the dashed curve in Figure 13.51 is almost indistinguishable 
from the numerically integrated solid curve. 

4. Lyapunov Exponents for Tw^o Planet Systems w^ith Mass Loss 

The discussion thus far has focused on single planet systems, whereas many solar systems 
contain multiple planets. In order to see how multiple planets affect orbital evolution during mass 
loss, we generalize the treatment to study systems consisting of two planets and a central star 
with decreasing mass. Such a 3-body configuration represents a crude model for our Solar System, 
where the motions of only the three most dominant objects (Jupiter, Saturn and the Sun) are 
considered. As a starting point, we fix the planetary masses and initial orbital elements (a, e) to 
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Fig. 6. — Value of the time/mass variable u = 1/m at the moment when the system becomes 
unbound, shown here as a function of the mass loss parameter 7, for varying values of the index /3. 
All of the curves correspond to circular starting orbits with rj = 1 (which corresponds to starting 
eccentricity e = 0). The curves correspond to values of /? = (top curve) to /3 = 3 (bottom curve). 
Dashed (black) lines correspond to the fitting function of equation (j79]) . 
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Fig. 7. — Value of the time scale ratio A evaluated at the moment when the system becomes 
unbound, shown here as a function of the mass loss parameter 7, for varying values of the index 
/3. All of the curves correspond to circular starting orbits r/ = 1 (which corresponds to starting 
eccentricity e = 0). The curves correspond to values of /3 = —0.5 (bottom curve) to 3.0 (top curve). 
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Fig. 8. — Value of the time scale ratio A evaluated at the moment when the system becomes 
unbound, shown here as a function of the mass loss index P which defines the time dependence of 
stellar mass loss. These values correspond to the limit of small mass loss parameter 7—^0 and the 
limiting case where the eccentricity of the starting orbit e = 0. The dashed curve, which is nearly 
identical to the solid curve, shows a simple fit to the function \f{l3), as described in the text. As 
shown in the text, A/ = 1 for the particular case 13 = 0. 
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those of Jupiter and Saturn, and set the initial stellar mass to Mq* = 1.0 Mq. We also restrict the 
orbits to a plane, thereby reducing the number of phase space variables from 18 to 12. To start, 
the mass loss function is taken to be an exponential model with index /3 = 1, although this law is 
generalized later. 

As long as the system suffers no close encounters, the orbit of an individual planet is similar 
to that described by the variable mass two-body problem. An example of the evolution of the 
osculating orbital elements (a, e) for our benchmark system (see above) is shown in Figure H] for 
a mass loss time scale of 10^ yr. As each planet orbits in an outward spiral, the semimajor axis 
increases approximately exponentially in time (in inverse proportion to the stellar mass). The 
eccentricity oscillates rapidly on orbital time scales, and more slowly on secular time scales (as the 
planets exchange angular momentum), but remains close to its starting value until stellar mass 
loss has taken place for a few e-folding times. The product of the semimajor axis and stellar mass 
(aM^) is approximately constant until a few e-folding times have elapsed. After a critical amount 
of mass is lost, the orbital elements a — )• oo and e — )• 1, and planets can become unbound. Notice 
that at this point, the stellar mass is only a few percent of the initial value; as a result, this scenario 
is rather artificial for stars like our Sun, which are only expected to lose about half of their initial 
masses. However, larger stars lose a greater fraction of their original masses. For example, a star 
with initial mass Mq* ~ 8Mq is expected to end its life as a white dwarf with roughly ~ 15% of 
its original mass (where the final mass fraction depends on the stellar metallicity). 

The evolution of the orbital elements can differ dramatically if the planets reach a small 
enough separation so that orbital crossings can occur. In this regime, chaos dominates and the 
orbital elements evolve in a less predictable manner. An example is shown in Figured] for the same 
parameters used in Figure [H but with the initial eccentricity of the inner planet (Jupiter) increased 
to e = 0.3 (a typical value for the current exoplanet sample). Although the system initially exhibits 
nearly periodic behavior, by the time t = t this stable evolution has been compromised. 

Studies of dynamical systems g enerally use the maximum Lyapunov exponent a s an indication 



of the level of chaos present (e.g., iLichtenberg &: Lieberman 



1983 



Strogatzj Il994l ) . If a system 



is chaotic, two nearby trajectories in phase space initially differing by a small amount Sq should 
diverge according to 



Sit) = 5oe 



At 



with 



A > 0. 



(81) 



The Lyapunov time is thus riy = 1/A. The long-t erm dynamical stability of the solar sy stem 
has been explored in the absence of stellar evolution (jBatvgin &: Laughlinll2008l : lLaskaiil2008l ). and 
current estimates of the L yapunov time for the solar system (while the mass of the Sun remains 
constant) are r ~ 5 Myr (ISussman &: Wisdomlll992l ). but this value decreases when stellar mass 
loss is introduced, as demonstrated below. 

Here we determine the Lyapunov times as a function of the mass loss time scale via numerical 
integrations. We define a "real" system along with a "shadow" system where the initial conditions 
of the shadow system differ by a small amount Sq. By integrating both systems simultaneously and 
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Fig. 9. — Osculating semimajor axis and orbital eccentricity for a pair of planets orbiting an initially 
solar-mass star with mass loss time scale r = 10^ years. Planets have masses, initial semimajor 
axis and eccentricities of Jupiter and Saturn. The orbital elements evolve in a roughly predictable 
manner, with the semimajor axes increasing smoothly and the eccentricities oscillating on secular 
time scales, but remaining relatively constant until the star has lost the majority of its initial mass. 



- 31 - 




Fig. 10. — Same as Figured! but with the initial eccentricity of Jupiter increased to eo = 0.3. This 
increase in eccentricity allows for orbital crossings and increases chaotic behavior. 
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monitoring the quantity 5{t), we can calculate the divergence of the neighboring trajectories and 
then estimate the maximum Lyapunov exponent. Since a three-body system restricted to a plane 
consists of 12 phase space variables, there is some choice in defining the quantity 6. For the sake 
of definiteness, we define 5 according to 



S = ^/{Xr-XsY + {yr-ys?, (82) 

where (x, y) are Cartesian coordinates for a planet's location, and where the subscripts r and s 
refer to the "real" and "shadow" trajectories respectively. Since chaotic systems display complicated 
behavior, the functions 5{t) will vary for effectively equivalent cases. As a result, for each system 
of interest, we run 1000 cases with both a "real" and a "shadow" system. To extract the Lyapunov 
exponent, we can either average together the 1000 runs to construct a single function 6{t) and use 
the result to find the exponent, or, we can find the exponent from each of the 1000 individual cases 
and then find the average exponent. Both schemes produce the same values; here we present results 
for the former case. Figure H] shows an example of the time evolution of 5{t) for different values 
of the mass loss time scale r. After an initial period of transient growth (roughly delimited by 
t < 0.3r), the divergence metric 6{t) increases exponentially with time and the Lyapunov exponent 
can be obtained by finding the slope of the line defined by In 5 = In 5q + At. 

For each value of the mass loss time scale r, the maximum Lyapunov exponent was calculated as 
outlined above. Since the maximum possible separation between the reference and shadow systems 
is finite (using the definition of 6 in equation [82]), the curves of growth eventually saturate. Thus, 
to extract the Lyapunov exponent A, we want to measure the curves of divergence after the initial 
interval of transient behavior but before saturation occurs. In most cases, A was calculated from 
the time-series data for times t/3 < t < r; this time interval is delimited by the vertical dashed 
lines in Figure HI An exception was made for the case of extremely slow mass loss, however, where 
r = 100 Myr. For this scenario, the time scale for mass loss is longer than the "natural" Lyapunov 
time of ~ 10 Myr (the value obtained without mass loss), and the curves of divergence saturate 
before t = 0.3r. In this case, A was calculated only for time series data with t < 10 Myr. Notice 
that the curves shown in Figure H] are not perfectly straight in the region between the dashed lines; 
this curvature introduces some uncertainty in the specification of the Lyapunov time scales. To 
estimate this uncertainty, we have calculated the Lyapunov time values for a wide range of choices 
for the time intervals. This procedure implies an uncertainty of a few percent. 

Figure H] shows our numerical values of the Lyapunov times riy as a function of the mass loss 
time r. We performed the analysis described above for each of the two planets in the system 
separately and averaged the results. The horizontal dotted line - included here for reference - 
corresponds to our numerically determined Lyapunov time for the model solar system in the absence 
of stellar mass loss. T his value is in relatively go od agreement with previous calculations for the 



complete Solar System (jSussman &: Wisdomlll992l ) , but differs slightly because our model considers 



only two of the four giant planets. Note that as r — )■ oo, the Lyapunov time approaches the dotted 
line, i.e., the value expected with no mass loss. The solid curve shows the calculated values of the 
Lyapunov time, whereas the dashed line indicates the least-squares fit (where the fit was taken over 
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Fig. 11. — Curves showing the divergence of two trajectories separated by a small distance (^o for 
different values of the mass loss time scale r. Black curve is r = 10^; blue is r = 10^; purple 
is r = 10^; and red is r = 10'' (yr). After a period of initial growth, the trajectories diverge 
exponentially, indicated by the linear shape of the latter portions of the graphs. The regions 
between the dashed lines were used to calculate the Lyapunov exponent. Note that the time 
variable has been scaled by the mass loss time scale r. 



- 34 - 



range of mass loss time scales r < 10 ). Notice that the slope of this line is close to unity. More 
specifically, we obtained 

riy ~ where p = 0.99. (83) 

Note that this fitted line cannot be meaningfully extrapolated below r = 10^ — 10^. In this regime, 
the mass loss time scale r becomes comparable to the orbital periods of the planets, and the 
dynamics of even single-planet systems becomes complicated (see the previous section). 

As a consistency check, we also explored other choices for the metric 5 that measures the 
difference between nearby trajectories. An especially compelling option is to use the semimajor 
axis a, because unlike the physical distance between the reference and shadow trajectories, there is 
no upper limit on this quantity. The previous calculation was thus repeated using 

S = \ar - Os|. (84) 

Our results are similar, however, which indicates that the Lyapunov times do not depend sensitively 
on the choice of 5. 

Next we would like to ensure that the (nearly) linear relation between the Lyapunov time and 
mass loss time is not an artifact of the exponential (/3 = 1) mass loss law that was chosen. Toward 
this end, we have explored two additional functional forms for the mass loss law: The first used 
vanishing mass loss index /? = (see equation [l9]), whereas the second used a constant mass loss 
rate with /3 = 2 (see equation [E]). For both of these mass loss functions, we obtained nearly the 
same power-law relation for the Lyapnunov time scale versus the mass loss time scale, i.e., T\y ~ r^, 
where p = 0.98 and p = 1.01 using equations ([19]) and (fT7|) . respectively. The results, shown in 
Figure HI are thus nearly identical, independent of the index (3 of the mass loss function. This 
finding suggests that this power-law trend is robust. 

Since the Lyapunov time scale is found to be comparable to the time scale for mass loss, we 
generally expect such solar systems to be only moderately influenced by chaos. In order for chaos to 
fully erase initial conditions for a dynamical system, nearby trajectories in phase space must diverge 
for several Lyapunov times. For example, in order for a starting uncertainty of 1° in position angle 
to grow into 360°, one needs ~ 6 Lyapunov times, or, about 6 mass loss time scales. For exponential 
mass loss, e.g., this time interval would result in the star losing 99.75% of its initial mass. Since 
most stars do not lose such a large percentage of their mass, the effects of chaos are not expected 
to completely erase the initial conditions of these systems. Nonetheless, chaos will partially erase 
the initial conditions; this trend will affect our ability to predict the phase of the orbit at the end of 
mass loss and will thus introduce uncertainty into predictions of the final (post-mass-loss) orbital 
elements (see Section 3.4). 
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Fig. 12. — Calculated Lyapunov times T\y versus mass loss time scales r for a star losing mass 
exponentially in time (mass loss index (3 = 1). The calculated quantities are shown (solid curve) 
along with a least-squares fit for r < 10^ yr (dashed line). As r — t- oo, the Lyapunov time 
approaches that of the solar system with constant stellar mass (riy ~ 5 Myr, as marked by the 
horizontal dotted line). 
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Fig. 13. — Lyapunov time scale t\j versus mass loss time scale r for systems where the star loses 
mass through three different decay laws. From top to bottom, the blue curve shows the results for 
/3 = 0; the black curve shows the results for exponential mass loss (/3 = 1); and the purple curve 
shows the results for constant mass loss rate (/3 = 2). For all three cases, the square symbols show 
the numerically calculated time scales and the dashed lines show a least-squares fit. These three 
examples show similar behavior, which indicates that the linear relationship between Lyapunov 
time scales and mass loss time scales is largely independent of the particular mass loss formula. 
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5. Applications 

To illustrate the efficacy of the results found in the previous sections, we consider two repre- 
sentative astronomical problems. For solar type stars, we find an effective outer edge of the solar 
system, i.e., the boundary between planetary bodies that remain bound after the epoch of mass 
loss and those that escape (Section 15. Ih . Next we consider planets that remain bound to white 
dwarfs, and find their final orbital elements (Section 15. 2p . 



5.1. Outer Boundary of the Solar System 



Suppose that a solar- type star, with initial mass Mq* = 1 M©, loses some portion of its mass 
over a time interval At = 1 Myr, so that the fraction ruf remains afterward. Since solar type stars 
lose most of their mass while they are on either the Red Giant Branch or Asymptotic Giant Branch 
(see Hurley et al. 2000), we expect /3 to be in the range 1-3. For a fixed time interval At and 
arbitrary index /3 > 1, the mass loss parameter 7 is given by 
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where qq is the initial semimajor axis of the orbit. For initially circular orbits, we can use the results 
of Section 3.5 to find the conditions for which planets become unbound. Using equation (j79p to 
define the critical value of 7, and equating the resulting value to the expression from equation ([85]) . 
we can solve for the critical semimajor axis Oc, such that orbits with larger values of a become 
unbound during the mass loss epoch. First we define the constant 
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The critical semimajor axis Oc is then given by 



a, = {Atf/^{GMo,)^/^A, 



m 



2/3 



(34,000AU)^/; 



m 



1 — m' 



(87) 



2/3 



1 



m 



13-1 
f 



We can thus find the critical semimajor axis ac for any given value of the index (3 and the remaining 
mass fraction mf. Note that the resulting value of Uc is a sensitive function of the index /3. To 
leading order, ac ~ (34,000 AU) m^^^"*^^^^^. If we use rnj = 1/2, the value expected for the Sun, 
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then the critical value of the semimajor axis Oc ~ 8500, 5350, and 3370 AU for indices /3 = 2, 
3, and 4 (these values for ac are in general agreement with the results of Veras & Wyatt 2012). 
For planets that escape, the corresponding velocities are small, in the range 0.3 ~ 0.5 km/s. For 
initially eccentric orbits, planets will become unbound for a wide range of starting semimajor axes, 
depending on the initial phase of the orbit (see Section [3]); however, this range is centered on the 
mean values found here. 

Note that in the limit /3 — t- oo, we formally obtain Oc — ?• for any value of the remaining 
mass fraction ruf < 1. However, this formulation of the problem breaks down before that limit is 
reached: For large values of /3, even though the time interval is fixed, the mass loss rate accelerates 
rapidly so that most of the mass is lost near the end of the time interval. Stellar mass is thus lost 
(effectively) through a step function in the limit of large /3. In this limit, the results of Section 3.3 
apply, and circular orbits become unbound (remain bound) for mass fraction ruf < 1/2 (m/ > 1/2). 

5.2. Orbital Elements for White Dwarf Planets 

Now consider a progenitor star with initial mass Mq* = 5 Mq, which evolves into a white 
dwarf with mass M^^ = 1 Mq (Hurley et al. 2000); the final mass fraction mj = 1/5 (uj = 5). For 
purposes of illustration, we assume that the mass is lost over a single epoch that can be described 
by a single value of the mass loss index /3, with a time scale At = 1 Myr. For orbits with starting 
semimajor axes a and eccentricities e, we would like to know the final orbital elements, after the 
epoch of mass loss. For orbits with starting semimajor axis a in the range 1 — 100 AU (closer 
planets are often accreted by the star), the value of 7 falls in the range 7 = 10~^ - 10-"^. The 
mass loss parameter 7 is thus small and nearly independent of the index /3 (see equation [85] ) . The 
time scale ratio A = ^u^~^^ f^^"^ . For the systems of interest, the largest 7 value is thus ~ 10~^, the 
largest value of u = 5, and the largest value of /3 = 3; the largest value of the time scale ratio is 
thus A ~ 0.063, so that A^ < 0.004. In the approximation scheme developed in Section 3.5, we have 
exact results when the integral J is small, where J = 0{X'^) < 0.004 (see also Appendix |A|) . 

The initial conditions for a planetary orbit include not only the semimajor axis and eccentricity 
(a, e), but also the phase of the orbit at t = 00 This latter quantity is specified by the initial value 
of the radial coordinate ^0 = /o; which lies in the range 1 — e<^o ^ 1 + e. With (,o{fo) determined, 
the integration constant E is given by equation (63). 

In the limit J — )• 0, the equation of motion shows that the function / oscillates back and 
forth between its turning points (given by equation [38]) while the system loses mass (analogous 
to the evolution depicted in Figure [33]) . If the duration of the mass loss epoch is specified exactly 
( equivalent ly, ifmj = l/uj is known exactly) , then we can determine the final value of the function 



full specification would also include three additional angles, e.g., the longitude of periastron, the longitude of 
the ascending node, and the inclination angle, but we can orient the coordinate system to eliminate this complication. 
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// = /(^/)) where Xf = uj. Let Nc denote the number of complete cycles that that function f{u) 
executes during the mass loss phase, where a cycle is defined as motion from one turning point to 
the other (a half orbit). In addition, the orbit must turn through a partial cycle from its starting 
value /o to the first turning point fj (where j = 1,2) and must turn through another partial cycle 
from the final turning point fk (where k = 1,2) to the final value ff. After integrating the equation 
of motion (|6Tll , we thus obtain 

_£;3/2 

[1 - m^] = {5N)o + N, + {5N) / , (88) 

where the integration constant E is given by equation (1370 and we have defined 

and 

where the it signs are chosen to keep the integrals positive. Note that equations (|88I- I9U|) completely 
specify the final value of the function ff. After solving these equations for /j, we can use equation 
(f66ll to find the final value <?/ of the orbital energy, and then use equation ([77|l and (iTHjl to find the 
semimajor axis and eccentricity of the orbit after mass loss has ended. 

Although the procedure described above is exact in principle (subject to the approximation 
J — )• 0), there exists a problem: In equation (f88]) . the quantities {6N)q and {6N) f are less than unity 
by definition, whereas Nc (and the left-hand side of the equation) is much larger, i.e., iVc ~ lO^-lO^ 
for the orbits considered here. To find the value // necessary to specify the phase of the final orbit, 
this approximate description for mass loss must be correct to better than 1 part in 10^ (10^) 
for orbits with starting a = 100 AU (1 AU). It is unlikely that the mass loss function for a real 
astronomical system obeys this simple model to such a high degree of fidelity. As result, even 
though we have an exact solution for the model equation, we cannot predict with certainty the 
final phase of the orbit for realistic systems. 

Given the uncertainty outlined above, post-mass-loss orbits can be described in terms of the 
expected values of the energy £f, semimajor axis aj, and eccentricity e/, as well as the possible 
ranges of values for the orbital elements (given the range of phases). For the sake of definiteness, 
we assume that the orbit starts at periastron (at the start of the mass loss epoch) and that the 
mass loss index /3 = 3. The orbit has starting angular momentum rj = 1 — e^. The integration 
constant E is then given by = 1 — 7^(1 — e)^ (from equation [36]) and the energy £j at the end 
of the mass loss epoch (from equation [75]) becomes 

= ^ {-1 + 390, 6257^ + (1 - 6)^7^ + 0{^^)} . (91) 
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Note that we can ignore the second 7^ term in the above expression. The expected final value of 
the (dimensionless) semimajor axis (from equation |77| ) is given by 

aj = 5 [1 -390,6257^]"^ , (92) 

and the corresponding expected value of the eccentricity (from equation [78]) is given by 

= + 390, 625r;7^ , (93) 

where we work to the same order of approximation as for a/ (recall that e is the starting, pre-mass- 
loss value of the eccentricity). Since 7 ~ 10~'^(ao/lAU)'^/^, the correction term is small: 390,625 
7^ ~ 4x 10~^(ao/lAU)'^ < 0.004 since oq < 100 AU. The mean value of final semimajor axis 
is thus about 5 times the starting value, as expected, and the leading order correction has been 
quantified. The square of the eccentricity increases by a similar amount. Because of the range of 
possible orbital phases at the end of mass loss, the orbital elements can differ from these mean 
values according to the relations 

Aa; _ A(e|) _ ^ , 125076. 

af e) £f ^ ^3/2 ' ^y^J 

where we have used equation (j76p and where el = + 7^(1 — e)^ is the effective eccentricity of the 
function f{u) during the epoch of mass loss. Thus, the total (relative) width of the distribution 
of possible final orbital elements is thus ~ 250076. Note that this range is often larger than the 
correction to the mean values. Consider a planet with starting semimajor axis = 100 AU and 
eccentricity e = 0.30. The mean value of the final semimajor axis is a/ ~ 502 AU, only 2 AU larger 
than the value suggested by the simple scaling law aM* w constant that is often used. However, the 
range of possible values about this mean is about Ao/ = ±19 AU. Similarly, the final eccentricity 
has mean value 6/ ~ 0.306, and the width of the range is about ±0.006. 



6. Conclusion 

6.1. Summary of Results 

This paper has reexamined the classic problem of the evolution of planetary orbits in the 
presence of stellar mass loss. Although this issue has been addressed in previous studies (see 
Section [T]) , we generalize existing work to include a new analytic formulation for time-dependent 
mass loss and to determine Lyapunov time scales for multiple planet systems. In particular, we 
consider a class of model equations where the mass loss index f3 is constant (see equation |10j). 
which allows for a wide range of time dependence for the mass loss rates and allows for a number 
of new results to be obtained analytically. Our main results can be summarized as follows: 

[1] Previous numerical studies show that planetary orbits often obey the approximate law 
am ~ constant, where a is the semimajor axis, and where this approximation holds as long as the 
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time scale for mass loss is significantly longer than the orbital period. By writing the equation of 
motion in the form given by ()22p and ()23p . we show analytically why this law holds (see equations 
[65] -[68]). In addition, the differential equation (jl2p for the energy E shows that the energy is a 
strictly increasing function of time, so that the semimajor axis (defined via a ~ l/l"?!) increases 
monotonically (whereas the orbital radius ^ oscillates in and out). 

[2] Previous literature often claims that the orbital eccentricity remains constant during the 
early phases of stellar mass loss. In contrast, this work shows that the eccentricity oscillates back 
and forth between well defined limits during the phase of mass loss, and that the amplitude of 
these oscillations grow with time (one example is shown in Figure [37^ . Moreover, the upper and 
lower limits of the eccentricity range can be calculated analytically using equations ([75]-[78|). Note 
that these oscillations in the eccentricity, while technically correct, result from assigning orbital 
elements (which describe ellipses) to orbital paths that are not elliptical. The actual orbit expands 
with time, and, for example, the outer turning point of the orbit increases monotonically (it does 
not oscillate). 

[3] In the limit of rapid mass loss, A — >• oo, we obtain analytic solutions that describe orbits 
for the entire class of mass loss functions (see Section [3. 3p . The condition for the planet becoming 
unbound is given by equation (I49p . For planets that remain bound, the new orbital elements are 
given by equation (f55]l . 

[4] Not all mass loss functions lead to planets becoming unbound (except, of course, in the 
extreme case where the stellar mass vanishes m — t- 0). The critical value of the mass loss index is 
/? = — 1, where systems with mass loss characterized by /3 < — 1 only lose planets in the m — ?■ limit. 
Note that systems with the transition value of the mass loss index /? = — 1 were first considered by 
Jeans (1924). 

[5] For the particular, intermediate value of the mass loss index /3 = 0, we can find analytic 
expressions for the function /(tt) and for the final values of the time scale ratio Aj when the planet 
becomes unbound (see Section [3^2]) . For arbitrary values of the mass loss index /? / 0, this approach 
can be generalized to find analytic expressions for f{u) and the orbital elements of the planet (see 
Section 13. 4p . The resulting expressions are approximate, correct to order ©(A^), and are thus 
accurate over most of the mass loss epoch. 

[6] One way to characterize the dynamics of these systems is through the parameter A / . We 
define A to be the ratio of dimensionless mass loss rate to the orbital frequency, and Aj is the value 
at the moment when the planet becomes unbound. For initially circular orbits, the parameter Ay- 
is always of order unity, and approaches a constant value in the limit of small mass loss rates 7; 
further, the value of this constant varies slowly with varying /? (see Figure [33]). For orbits starting 
with nonzero eccentricity, however, the parameter Aj can depart substantially from unity and varies 
significantly with /3 (compare Figures 13.21 and l3.2p . 

[7] For multiple planet systems, we find that the Lyapunov times decrease in the presence of 
stellar mass loss, so that chaos should play a larger role in planetary dynamics as stars leave the 
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main sequence. In fact, the Lyapunov time scale is proportional to (but somewhat shorter than) 
the mass loss time scale over a range of conditions (see Figures H] and H]) . For a typical mass loss 
time scale of r ~ 10^ yr, the Lyapunov time is ~ 2 — 3 x 10^ yr. Three different forms for the stellar 
mass loss function have been considered, all yielding similar results, which suggests that this trend 
is robust. 

6.2. Discussion 

In spite of its apparent simplicity, the classic problem of planetary orbits with stellar mass loss 
is dynamically rich. Here we discuss two issues that have been highlighted by this present work: 

The use of osculating orbital elements is a standard way to describe planetary motion. In this 
scheme, the planet is assigned the Keplerian orbital elements that it would have if it were moving 
within a purely Keplerian potential. In systems with stellar mass loss, however, this approach 
might be more misleading than illuminating (see also Radzievskii & Gel'Fgat 1957, Hadjidemetriou 
1963). Consider, for example, the osculating eccentricity of the orbit. As the star loses mass, 
the eccentricity oscillates (Figure 13. 4p with an amplitude that grows with time. This oscillating 
eccentricity would seem to imply motion along an elliptical path, where the shape of the ellipse 
cycles back and forth between being more elongated and more round. However, this "oscillation" of 
the eccentricity is an artifact of assigning a Keplerian orbital element (e) to orbital motion that is 
not Keplerian. During the mass loss epoch, the orbit has inner and outer turning points, analogous 
to those of an ellipse. But the turning points of the actual orbit do not oscillate - the orbit smoothly 
spirals outward (e.g., see Figure 2 of Veras et al. 2011). 

An alternate description of the dynamics can be constructed: Before the onset of mass loss, the 
orbit is a Keplerian ellipse and can be described by the usual orbital elements (a, e), which provide 
the initial conditions for the next stage (along with the phase of the orbit at t = 0). During the 
mass loss epoch, the orbit is not an ellipse, but the scaled radial function / = = £,/u follows a 
nearly Keplerian trajectory. For the particular mass loss function with index /3 = 0, this analogy is 
exact. In this case, the function / obeys the same equation of motion (|35]) as the radial coordinate 
^ in the Kepler problem (albeit with a nonconventional time variable where dt ~ du/v?). The 
scaled function / has turning points (equation [38]), an effective semimajor axis a* = 1/E (see 
equation [39]), and an effective eccentricity e* (equation [10]). In the general case where /3 7^ 0, 
the scaled function / executes nearly Keplerian motion: Here the equation of motion for / has a 
Keplerian form up to a correction of order 0{J) = 0{}?), where A ^ 1 for most of the evolution 
(until shortly before the planet becomes unbound) . Note that the effective orbital elements (a* , e* ) 
characterizing the function / are constant, to O(A^), while the star loses mass. After the epoch 
of mass loss, the planet (if it remains bound) once again enters into a Keplerian orbit, now with 
elements (aj,ej). The problem can thus be described in terms of three sets of orbital elements: 
before (a, e), during (a*, e*), and after (a/, e/) mass loss. But the orbital elements during mass loss 
(a^,,e*) correspond to orbits of the scaled function / (the orbit in physical space is not Keplerian). 
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Another interesting complication arises: These orbits display a type of sensitive dependence 
on initial conditions even in the absence of chaos. For planets that remain bound after the mass 
loss epoch, the final orbital elements (a/,e/) depend on the phase of the orbit, at both the start 
and the end of the mass loss epoch. However, the number of orbital cycles during mass loss is large, 
A'^c ~ 1/7 ^ 1, and a precision of one part in A^ is necessary to specify the final orbital phase (given 
the initial phase). The orbit can be calculated to sufficient accuracy to specify the final orbital 
elements provided that the starting state, the duration of stellar mass loss, and the form of the 
mass loss function are all known well enough. In practice, however, real astronomical systems will 
not follow these particular forms to such high accuracy, so that the final orbital elements cannot be 
predicted with certainty. Instead, the expectation value of the orbital elements can be predicted, 
along with the range of possible variation about their mean values (see also Section [O]) . 



6.3. Future Work 

There are many opportunities for future work. Analytic studies can be taken into two di- 
rections. First, the formulation developed here can be applied to a wide range of astronomical 
systems, including predictions of orbital elements for planets that remain bound after stellar mass 
loss and predictions of the conditions required for bodies to become unbound. On the other hand, 
the analytic treatment can be developed further to include more general mass loss functions, mul- 
tiple phases of mass loss, and higher order approximations for the conditions under which planets 
become unbound. 

For the calculations of the Lyapunov time scales, this paper has focused on analogs of our own 
solar system, and has considered only the motion of Jupiter, Saturn, and the Sun, since these are the 
most gravitationally dominant bodies. However, future calculations should also include additional 
planets and a wider range of starting orbital elements. In addition, this paper has considered 
relatively short integrations, spanning at most only 10 — 100 Myr. Although stellar mass loss is 
not expected to continue for longer than 100 Myr, the increased semimajor axes of the planets will 
change the overall dynamics and the decreased stellar mass (relative to the planets) could lead to 
an increase in dynamical instabilities. As a result, longer integrations should be performed, using 
the lower (constant) stellar mass and the increased semimajor axes of the planets as input. The 
calculations of this paper are limited to coplanar planeteary systems; future work should explore 
the effects of different inclination angles. Finally, given the diversity exhibited in the observed 
sample of extrasolar planets, different planetary configurations should also be considered. These 
types of calculations will help us understand the long term fate of planetary systems in general and 
will help direct future observations. 
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A. Bounds on the J Integral 

In the limit \J\ <C 1, we have a completely analytic description of the dynamics. It is thus 
useful to place bounds on the integral J, which can be done as follows: First write 

J = -/'^I3I where ^^^/^ x^ffxdx. (Al) 
Note that since the function / is of order unity, the order of the integral J is given by 

J = (7^x2) = O (A2) . (A2) 

Next we integrate by parts to obtain 

I = x^f-fl-2j\fdx. (A3) 
The second integral can be written 

2 xfdx = 2{f) xdx = if) (x^ - 1) , (A4) 

where we have invoked the mean value theorem. The function / varies between turning points so 
that /i < / < /2- We thus have bounds 

I < Afl - fl) + fi - fo < Afi - fi) , (A5) 

and 

I > x\fl - /I) + /I - > Afi - /I) ■ (A6) 
As a result, we have the bound 

|/| < x^ifl - ff) = x\h + h){h - h) = xHale. . (A7) 

We thus obtain the desired bound on J, i.e., 

\J\ < M^x\f2 + /i)(/2 - /i) = l^xH^ale, . (A8) 

In order to evaluate this bound, we need expressions for the turning points /i and /2, or, equiva- 
lently, the effective semimajor axis a* and eccentricity e*. As derived in the text, we have approx- 
imations for these quantities, where these expressions are exact in the limit J — > 0. 
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